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We analyze the dynamics of a driven, damped pendulum as used in clocks. We derive equations 
for the amplitude A and the phase if of the oscillation, on time scales longer than the pendulum 
period. The equations are first order ODEs and permit fast simulations of the joint effects of 
circular and escapement errors and friction for long times (1 year real time in a few minutes). The 
equations contain two averages of the driving torque over a period, so that the results are not very 
sensitive to the ‘fine structure’ of the driving. We adopt a constant-torque escapement and study 
(1) the stationary pendulum rate ip as a, function of driving torque and friction, (2) the reaction of 
the pendulum to a sudden change in the driving torque, and (3) to stationary noisy driving. The 
equations for A and ip are shown to describe the pendulum dynamics quite well on time scales of a 
period and much longer. The emphasis is on a clear exposition of the physics. 


I. INTRODUCTION 

Suppose you wish to study the long-term time keeping 
properties of that longcase clock that for years on end 
has been peacefully eating away the seconds of your life. 
How would you do that? Making regular measurements 
is one thing, and numerical simulations is another. In the 
latter case you are in for a surprise: you have to solve the 
equation of a driven, damped harmonic oscillator in steps 
much smaller than the period of the pendulum, allowing 
for the slow changes (in driving, friction, temperature,..). 
And although each next oscillation is virtually identical 
to the previous one, you must follow them all in great 
detail to get any accuracy in the final result. This leads 
to long simulations, even in these days of fast desktop 
computers. Here we develop and validate a much faster 
method, that allows making time steps of the size of the 
pendulum period. 

There is an extensive literature on the various types of 
pendulum and the effect of perturbations on their mo¬ 
tion. Rawlingsi and Woodward^ expound the scientific 
principles underlying clocks and pendulums, while Baker 
and Blackburn^ review the diversity of pendulums occur¬ 
ring in nature. A useful (but incomplete) compilation of 
the literature has been given by Gauld4 

The oldest known pendulum timing error is the circu¬ 
lar error, i.e. the fact that the period of a free pendulum 
depends weakly on its amplitude. Amplitude variations 
induce therefore small timing errors. Huygens showed in 
1657 how these may be eliminated by the use of cycloidal 
cheeksHenceforth we shall often just speak of errors in¬ 
stead of timing errors. The escapement,^ i.e. the mech¬ 
anism that transfers energy to the pendulum to keep it 
swinging, is an unavoidable source of errors. These have 
been computed by Airy in 1830 in a seminal paper that 
preludes all later developments on this topic ^ 

Kesteven^ computed the period change of a pendulum 
driven by a dead-beat or a gravity escapement, while 
Nelson and Olssoii^ did so for air drag and buoyancy. 
Denny^ shows how a simple escapement makes a pendu¬ 
lum settle in a stable limit cycle, and presents a useful 


summary of Airy’s method. Friction is often modelled as 
a force proportional to the velocity, but Andronov et al^ 
treat several other types of friction. 

Advanced models recognize the fact that a pendulum 
with its escapement is a mechanical system with several 
degrees of freedom that may be analyzed using so-called 
impulsive differential equations, see e.g. Moline et 
and Moon and Stiefeld^ On the practical side, Matthysi^ 
gives detailed suggestions for the design and construc¬ 
tion of a precision pendulum, with much attention to the 
properties and choice of materials. 

The basis for a study of a pendulum’s time keeping 
properties is the equation of motion. It incorporates 
all effects but requires a time step much smaller than 
the period and contains therefore too much information. 
Here we derive equations for the pendulum dynamics on 
a meso-time scale, i.e. the time scale of a period. This 
work amounts to a new, more general formulation of es¬ 
capement theory. It permits much faster simulations, 
while all timing error sources operate in concert. 

The technique is explained in Secs. [Ill and imi We in¬ 
troduce a simple model escapement in Sec. mm and 
analyse the stationary operation of a pendulum driven 
by this escapement in Sec. El In Secs. El and EH we 
apply the theory to study the reaction of the pendulum 
to a sudden parameter change and to stationary noisy 
driving. The second novel aspect of our work is that we 
validate our approach in detail by comparing the results 
with data computed directly from the equation of mo¬ 
tion. We discuss our results in Sec. IVHI 

II. EQUATIONS FOR THE AMPLITUDE AND 
PHASE 

The pendulum oscillates in a plane, making an angle 
ip{t) with the vertical. The equation of motion for (p is 

(fi + 2e(p + sintp = m{(p) . (1) 

Here ip is measured in radians, and e and m are the 
friction coefficient and the driving torque (divided by the 
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FIG. 1: The motion of the pendulum is represented by a 
point with radius vector A rotating with angular speed ipj 
projected on the vertical axis, so that ip = Asin^. In this 
way we achieve that (p = 0 when xfi is an integer multiple of tt. 
Both A and xli exhibit slow variations (slower than indicated), 
resulting in gradual changes in the swing amplitude and the 
period of the pendulum. 

moment of inertia of the pendulum). The dot stands for 
the time-derivative: ' = d/dt, "= d^/dt^, etc. Time has 
been scaled so that the period of unperturbed, small- 
amplitude oscillations is 27r (e = m = 0; sin<p ~ if). We 
shall refer to this as the nominal period. 

The escapement^ is designed to deliver a torque m that 
depends only on the angle f. It has no intrinsic fre¬ 
quency of its own, but just delivers the torque required 
by the value of f, irrespective of time. So, m = rn(ip) 
and Eq. m is autonomous. In actual fact m may have 
a weak dependence on ip as well. We shall not consider 
the more general case m = m{ip, ip), and we shall take e 
constant. But the theory developed below can handle e 
and m with arbitrary functional dependence on ip and ip. 

It is straightforward to solve Eq. o numerically, but 
the result contains too much information for our pur¬ 
poses. The time keeping is determined by the times of 
successive = 0 crossings, for which it is sufficient to 
know ip{t) with a time resolution of about a period. In 
other words, we are interested in slow changes of the 
amplitude and phase. The problem of gradual secular 
changes in quasi-periodic systems has come up in physics 
time and again. One of the earliest applications was the 
problem of secular perturbations in planetary orbitsJ^ 
Here we shall follow the theory developed by Krylov, Bo- 
goliubov and Mitropolski (KBM) The idea is to take 

if = Asin^ ; p = \ cosip , (2) 

where A is the maximum of p (referred to as ‘the ampli¬ 
tude’), and Ip the phase of the pendulum, see Fig. [T] 
We may regard ip as the time measured by the pen¬ 
dulum (‘pendulum time’), while t is a stable reference 
time (‘laboratory time’). One might think that the 
first relation p = Xsinip suffices, for we may compute 
ip = X sin Ip + Xip cos Ip and then ip =■■■, insert all that in 
Eq. ([T]) and make approximations using A, e I. How¬ 
ever, that path is fraught with ambiguities due to the 


fact that p = X sin ip does not define A and ip uniquely. 
We do need a second relation. 


A. The method of Krylov, Bogoliubov and 
Mitropolski (KBM) 


The next i 

step is to solve Eq. 

(ID) for A and ip, 


A = 


■ lP 

= arctan((^/(/:) , 

(3) 

and to compute the derivatives: 



A = 

pp pip 
(<^2 + (^2)1/2 

' -( 1 

II 

{p + p)p , 

(4) 

iP = 

p^ - pip 
p"^ + p^ 

= 1 - 

-^i‘P + p)p . 

(5) 


The derivative ip of the phase is the rate at which the 
pendulum ticks. And ^A^ is the total (kinetic -I- po¬ 
tential) energy of the pendulum. Both right hand sides 
contain p + p, which the equation of motion (HD allows 
us to write as: 

p + p = —2ep + (p — sin p) + m{p) . (6) 

On the right we recognize the terms responsible for the 
effects of friction, circular error and driving, respectively. 
We insert Eq. ([61 in Eqs. m and ©, and eliminate the 
remaining p and p with the help of Eq. 

A = —2eAcos^ V'+ {•^sini/j — sin(Asin^)} cos'f/' 

-bmcos^, (7) 

Ip — 1 = 2ecos'!/'sin'!/;-|-{A“^ sin(Asin'0) — sin^/:} sin-^ 

— (m/X) sinip . (8) 

Equations m and m are still equivalent to Eq. (HD. 
KBM then average over one period to remove the fast 
time scales. This amounts to integrating the equa¬ 
tions over Ip over the interval [^o “ V’o + ^] centered 
at an arbitrary ipQ, and to assume that e, X and ip 
but not m, are constant over one period. Afterwards, 
the index 0 on '0 is dropped. In doing so, the terms 
{Asin^/’ — sin(A sin'0)} cosip in Eq. ([7D and 2e cosip sin ip 
in Eq. m vanish because they are antisymmetric in ip. 
Since {2tt)~^ cos^ ip dip = ^ we obtain: 

A = —eX + j) m{p) cos Ip dip + Oi(e^) , (9) 

ip — l = -—-- (f m{p) sin Ip dip + 02 (e^) ,(10) 

16 A / 

where we use the notation f = (27r)~^ fZ.„.. We demon¬ 
strate in appendix lAl how the circular error term —A^/16 
is extracted from Eq. 
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In reality e, A and "0 are not constant over a period. 
Bogoliubov and Mitropolskii^ show that this produces 
correction terms in Eqs. and OT having the form of 
a power series in £. These higher order terms play no 
role here, except in Sec. EH where we believe to observe 
the effect of the lowest order terms oc £^. 

Briefly, what if e and m depend on ip and pi Take 
for example £. Then set £ = e{p,ip) = £(A sin-(/i, A cos (/?) 
in Eqs. © and ©. The average over a period will now 
generate a friction-like term in the rate equation OT, 
etc. 


III. THE PHYSICS OF EQS. dH) AND (ITUD 

Equations © and (uni) are first order ODEs well- 
suited for a study of the long-term dynamics of the pen¬ 
dulum. Equation © determines the amplitude A and 
the total energy A^/2 of the pendulum. The amplitude 
changes only slowly because the friction time scale 
is long, l/(27r£) nominal periods, and because the es¬ 
capement needs to deliver many small kicks to alter the 
energy ^A^. The damping is compensated by the driving 
term f mcosipd'tp. A back-of-the-envelope derivation of 
Eqs. © and (fTUl) is given in appendix [BJ 


A. The rate equation (I10|i 

Equation m for the rate features the circular error 
term — A^/16, due to the difference between p and sin(/j 
in Eq. ©. It has a minus sign because the circular error 
increases the period, so it must reduce the rate. Friction 
does not appear, and has therefore only an indirect effect 
on the rate through its influence on A. 

The term § m sin if) dtp is called the escapement 

error. It specifies the relative rate change induced by a 
driving torque m, as computed already by Airy^. The 
history of horology is to a large extent a quest for tech¬ 
nical solutions to make this term and the circular error 
as small as possible. 

An important feature of both equations is that the 
driving torque m appears only in the form of two aver¬ 
ages over a period. This is a consequence of the elim¬ 
ination of fast time scales and demonstrates that much 
of the fine structure in the driving torque averages out 
and is irrelevant for the pendulum’s long term dynam¬ 
ics. Only two numbers matter, namely the even (cosine) 
and odd (sine) moment, and all pendulums with escape¬ 
ments having the same j> m cos ip dip and j> m sin ip dip are 
in principle equivalent time keepers. 


B. A model escapement 

For our numerical work we introduce a simple escape¬ 
ment, see Fig. [5] and Woodwardi^. A constant driving 
torque m = p > 0 switching to m = —p after half a 



FIG. 2: The model escapement delivers a constant torque 
m = p, switching to m = —p after half a period (— — —). 
The switch is at a fixed pendulum angle po in the first quarter 
of the swing, and —po in the third. The phase ip is indicated 
at the bottom. 


period, and so on. The switch is at a fixed (mechanically 
determined) pendulum angle po- The phase ipo at that 
moment depends on the amplitude according to Eq. ©: 

y)o = Asin'(/)o ; 0 < V'o < ’’■/2 . (H) 


The swing amplitude A must be larger than po , otherwise 
the driving torque cannot switch sign, and the pendulum 
will stall. Since po = Asin^o has two solutions, ipo and 
TT — ipo, we impose 0 < ipo < tt/2 because escapements 
generally switch to a torque counteracting the motion 
before the pendulum reaches maximum amplitude. 

Actual escapements deliver a more complicated 
TO^ iiOiii but as argued above, this fine structure of m 
is expected to be relatively unimportant for the perfor¬ 
mance. So we hope to catch the main features of this and 
similar escapements with the help of this simple model. 
We compute the two averages required in Eqs. © and 

(US: 


a = 


m cos'^dV’ 


( 12 ) 


1 

27r 



2-K+'lpo 

TT-l-l/'O 


p COS Ip dip 


. I 

= — sin-^o 

TT 


‘2ppo 

ttA 


( 13 ) 


Since the integrand is a periodic function of ip we may 
start the integration where we like, as long as we integrate 
over a full period. A similar computation yields 


b = 


m sin ip dip 


2p 


cos "00 ■ 


( 14 ) 


Equations © and OT for the amplitude and the phase 
may now be written as: 


dt 


-2eX^ 


dppo 

TT 




( 15 ) 


Ip-l = ^ COSlpo + 02 (£^) , ( 16 ) 

lo TTA 
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where Eq. (1161) may be further reduced with cosijjo = 
[1 — as follows from Eq. (fTTl) . Due to a dif¬ 

ferent notation the escapement error does not look like 
any of the expressions in Ch. 7 of Rawlingai. But we 
have verified that we can reproduce these formulae start¬ 
ing from our expressions [2/j,/(7rA)] cosipo in Eq. (fTHll or 
— (1/A) j> TOsini/'dV' in Eq. (ITUl) . So the difference is only 
apparent. In appendix E] we analyse the action of an 
escapement that delivers a single impulse in each period. 


C. Analytic solution 

Before we start using Eqs. m and (HU), we point out 
that they often possess a simple analytic solution, which 
may be useful in theoretical investigations. For example, 
if e is constant and the transient responses are over, the 
system is in a quasi-stationary state and the solution of 
Eq. ((15]) is: 

^2 _ f exp(—2er)/r(t — t) dr . (17) 

Jo 

This permits us to compute A(t) for any time-dependent 
driving torque And once A is known, we may com¬ 
pute Ip from Eq. (HU. For example, let there be a jump 
in /r at t = 0 and fj, = for t < 0 and p, = /ii for t > 0. 
Inserting that in Eq. HID produces 

\ - (mi - Mo) exp(-2et)} t > 0 . 

( 18 ) 

This result will be used later to explain the numerical 
work in Sec. 0 Just replace t by t — to everywhere in 
Eq. (fT51) if the jump is at to- 


IV. STATIONARY OPERATION 

The asymptotic theory developed above enables us to 
address many questions that would be difficult to answer 
if we had only the equation of motion (HD. We begin with 
the case that e, /i and ifo are constant. The pendulum 
will then perform stationary oscillations after some time. 
We study how the stationary amplitude A and rate tp 
depend on the parameters. The amplitude follows by 
setting dX^/dt = 0 in Eq. (ITSl) : 

X = . (19) 

\Tr£(poJ 

and using this in Eq. HU we obtain the stationary rate: 



e 


FIG. 3: Contour plot of — 1 in sec/day as a function of e 
and fi, for ifio = 0.03 = 1.7°. This is the rate in excess of the 
nominal rate of a pendulum driven by the model escapement. 
Lines of constant swing amplitude are dotted, A = ipo and 
higher up 2(fio and 3(po. For example, if e ~ 5 x 10~® and 
H ~ 10“®, the pendulum loses about 10 s per day with respect 
to the nominal rate, and the swing amplitude is about 2ipo. 
Further details in main text. 


the rate. The second term is the escapement error, which 
is manifestly positive. It follows that the driving nec¬ 
essary to keep the pendulum going, always forces it to 
run a little faster as well. This confirms one’s intuition, 
and the physical reason is, loosely speaking, that during 
most of the time the escapement is pushing the pendu¬ 
lum ever so gently forward. These statements hold for 
the model escapement and comparable types, e.g. the 
recoil escapement^ but not for all, see appendix [Cl 
Fig. [3] shows the excess rate on top of the nominal rate 
^p = 1 in sec/day due to circular and driving error to¬ 
gether. The £-fi plane is divided in three regions: 

1. A losing rate, to the left and up. In this region the am¬ 
plitude A is large because the friction is relatively small. 
As a result, the circular error is larger than the escape¬ 
ment error, and the pendulum rate is sub-nominal. 

2. A gaining rate, to the right and up: here the effect of 
the driving torque dominates (escapement error). These 
regions are separated by the dashed contour where circu¬ 
lar and escapement error compensate each other. Pendu¬ 
lums operating on this dashed line run at nominal rate. 

3. In the third region, below the line A = <po, the driv¬ 
ing is too weak to keep the swing amplitude A above 
the value <po, and the pendulum stalls. The minimum 
required torque p, is ^7r£(po- 


A. Losing or gaining? 


iP-1 = --— -b e 
Sne 


2fj, 

TTEipo 



( 20 ) 


The first term on the right is the circular error (no longer 
readily recognizable as such). It is negative and reduces 


Here comes a perennial question: will a clock lose or 
gain when a parameter is changed? Take for example 
the driving torque. The discussion usually proceeds as 
follows: p. t implies a larger amplitude A, which means a 
larger circular error, i.e. ip j,. But sometimes we observe 
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the opposite. How come? On closer look, the answer de¬ 
pends on the pendulum’s initial operating point in Fig. [31 
If that is towards the left where the circular error dom¬ 
inates, then /r makes the first term (circular error) in 
Eq. ([20jl go down and tp |. But if we start more to the 
right where the escapement error dominates, the second 
term in Eq. (1201) (escapement error) increases. So tp 
when we start on the left in Fig. [31 but ip t when we 
begin on the right. There is no simple answer. 

Variations in friction may be dealt with in the same 
manner. A remarkable feature is that the pendulum may 
run faster when friction is increased. Take for example /r 
fixed at 2 x 10“® in Fig. [31 then the rate goes up with e 
until e ^ 2.5 x 10“^ and decreases only thereafter. 


B. Accuracy 

The accuracy of the time keeping is a subject with 
many different aspects, of which we shall discuss only one. 
Experience has shown that the accuracy of a pendulum 
is largely set by the friction coefficient e. The driving 
plays a minor role. Why is that so? The answer comes 
in two steps. 

We focus on long time scales £~^)- In that case 
we may ignore the time derivative in Eq. and obtain 
the following order-of-magnitude estimate: eA ^ m, or 
mjX ^ e, which we use to estimate tp — 1 with Eq. (ITOl) : 

tp — 1 ^ m/X s . (21) 

Here we assume that the pendulum operates in region 
2 of Fig. [31 where the circular error is negligible. Thus 
we have found that ip — 1 ^ e, and Eq. (EOl) is basically 
telling the same story. The point is that the driving 
torque no longer appears in Eq. (1311) because it could be 
approximately eliminated by requiring quasi-stationary 
operation. 

The second part of the answer is that as long as the 
parameters are constant, nothing changes. The rate is ex¬ 
actly constant and the pendulum is a perfect time keeper, 
even though it may not run at nominal rate. But param¬ 
eters are never constant and change on a variety of time 
scales. Experience has shown that changes of several tens 
of percent in the parameters and more are not uncom¬ 
mon. And that makes s a reasonable upper limit for the 
relative precision of the pendulum rate. 


C. Optimal driving torque 

From Eqs. m and (nni) the least disturbing driving 
torque is seen to be an impulse in a narrow phase interval 
6ip 1 centered at ^ = kn, when the pendulum passes 
through its equilibrium position ip = 0. The escapement 
error is then X~^ f msmipdip ~ {m/X){dip)'^ ~ e6ip. 
Here we used that m6ip ~ eA, as follows from Eq. ([^. 
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FIG. 4: Numerically determined pendulum angle p. On 
the horizontal axis time in nominal periods. Parameters: 
e = 2.5 X 10~^, ipo = 0.03, timestep dt = 0.03. Initial con¬ 
ditions: :/p(0) = 0, </5(0) = (2po</5o/7r£)^^^ (stationary swing). 
hXtl2'K = 150 the driving torque switches from /tq = 2 x 10“® 
to /ii = 4 X 10“®. There is no visible reaction. 

We conclude that it should be possible to beat the pre¬ 
cision estimate (EU by a factor 5ip with the help of 
a specially designed escapement. All this is well known - 
we merely illustrate that the present formulation which 
uses concepts that are close to the instrument, allows a 
more direct answer to such questions than Eq. alone. 

V. TRANSIENT EFFECTS 

We validate the asymptotic theory developed above 
with two case studies. The first is the reaction of the pen¬ 
dulum to a jump in the driving torque. We solve Eq. o 
numerically, parameters as in Fig. [4] for details see ap¬ 
pendix [D] The pendulum is set off in a stationary swing, 
and the value of 1/3(0) follows from Eq. ([3]): <^(0) = A and 
A may be had from Eq. (1131) . There are 2tt/ dt ^ 209 
time steps in one nominal period. A similar numerical 
experiment has been carried out by WoodwardjI^ 

The solution. Fig. 01 is not very telling: a series of 
identical oscillations, and there is no visible reaction in 
the vicinity of the jump in p, at t/27r = 150. To refine 
the analysis we determine the times U that tppt) crosses 
zero (from — to -f). Numerical details are given in ap¬ 
pendix [D] Since ip increases by 27r at each next zero i of 
ip (by definition), we may determine the time (i.e. phase) 
Ip and the rate ip measured by the pendulum on a non- 
equidistant grid {U}, as follows: 

ipi = 2-Ki ; ipi = 2-K[ti - . (22) 

This procedure can be easily implemented in the numer¬ 
ical model (and in a pendulum experiment). The mea¬ 
sured rates may then be confronted with theory, Eq. (El]). 

The reaction of the pendulum. Fig. [SJ top panel, is 
now clearly visible. The initial operating point of the 
pendulum is A in Fig. [31 where it gains about 1.1 x 10““^ 
or ^ 10 sec/day. The rate changes instantaneously, to¬ 
gether with the jump in the driving, and then tapers off 
to a new stationary value ^ 2 x 10“'* or ~ 17 s/day in 
the new operating point B. The time scale is set by A, 
Fig. [^ bottom panel, and A changes only slowly because 
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FIG. 5: Reaction of the pendulum to a jump in the driving 
torque (parameters as in Fig. |4]). On the horizontal axis time 
in nominal periods. At t/2TT — 150 the driving torque switches 
from /ro = 2 X 10~® to /ri = 4 x 10~®. Top panel: excess rate 
of the pendulum (in seconds per sec), as ‘observed’, that is, 
computed from dehnition (EH). The noise on the rate curve 
is numerical and can be reduced by taking a smaller timestep 
dt. Diamonds o are values inferred from theory, Eq. (HI. 
Bottom panel: amplitude A computed from Eq. m and the 
simulation data; Eq. (El) yields the same result. 


it takes many small pushes of the escapement to change 
the energy A^/2. Once A is known, we compute ip — l with 
Eq. (fTC)) . as a verification. A few values are displayed as 
o in Fig. [S] 

The magnitude of the rate jump follows directly from 
Eq. (fTCl) : since A changes slowly, and costpo depends only 
on A, the rate difference just before and after the jump is 
dtp = 2{6fi) costpo /{nX). This allows us to compute that 
tp — 1 = 3.2 X 10“^ after the jump, as observed. The jump 
in the rate is tiny but well measurable. It is an arresting 
experience for students experimenting with a pendulum 
in a practical physics course. The unaided eye sees no 
change when the driving is altered, and yet it happens, 
just because the escapement pushes a little harder. 

Finally, we compute the time shift the pendulum in¬ 
curs. Locating the jump for simplicity at t = 0, the phase 
shift is Atp = J^{tp — tpr) dt, where tp is the rate given by 
Eq. (fTBl) using Eq. (IT^ for A, and tp,- the constant rate of 
a reference pendulum operating in point B, i.e. the area 
under the curve tp and above the tangent to tp in t = oo. 
We find Atp = 0.295, or {Atp/2'K) • 2s ^ 0.1s for aim 
pendulum. 

We have also looked at the pendulum’s reaction to sud¬ 
den changes in e. The main difference is that the rate no 


longer makes a sudden jerk but adapts slowly to its fi¬ 
nal value, like A does. In conclusion we can say that the 
asymptotic theory seems to work well, also for a rapid 
change as we encountered here. 


VI. LONG-TERM BEHAVIOR 

The second case study is a pendulum with a variable 
driving torque, running at nominal rate in point C in 
Fig. [3] in the absence of variability. The co-ordinates of 
C are /tq = 1.5 x 10“^ and e = 1.167245 x 10“'* (obtained 
by solving Eq. (fT31) for e with tp = I and /tq = 1-5 x 10“®). 

The reason for selecting this model is not because we 
think it is particularly realistic, but because we need a 
‘rough’ test bed to see how well the data generated by the 
‘numerical clock’, i.e. Eqs. (HD and (l23|) , are recovered by 
the theoretical model, i.e. Eqs. m and m together 
with (1331) . 

The simulation advances Eq. dm for 10® timesteps 
{dt = 0.03), comprising ~ 4.8 x 10® nominal periods, 
or about 110 days for aim pendulum. The initial con¬ 
ditions are as in Fig. |3] (stationary swing). During the 
simulation we determine the —h zero crossing times {t^} 
of ip and the pendulum time tpp^i at these moments, ba¬ 
sically by counting sine waves, see definition dUD. We 
also integrate Eqs. ED and (fT6l) on this {U} grid, to 
determine the model pendulum time 'i/’m.i- 

The driving torque fj, has a random component: 

p. = fio{l + 6n{t)/p.o} , (23) 

where 6fi{t) is a normally-distributed random variable 
with zero mean and variance SpP^ that is renewed at 
every —I- zero crossing of ip. The coherence or correlation 
time of Sp.{t) is therefore one nominal period. 

In spite of the location of C in Fig. |31 the pendulum 
does not run at nominal rate. In the absence of variabil¬ 
ity it gains ^ 0.005/5 x 10“® = 10“®, or 0.3s/yr, see 
Fig. [51 top panel. This is ten times better than the per¬ 
formance of a Shortt clock, see Ref. [l|, p. 126. The effect 
is therefore extremely small, but we discuss it because we 
are after long-term effects here. Otherwise, the pendulum 
is seen to be almost impervious to the large variability in 
the driving, with an r.m.s. gain of ^ 0.04/4 x 10® ~ 10“®, 
or ^ 0.3s/yr (same number by coincidence). 

Why is it so insensitive to the fluctuations? The deci¬ 
sive factor is the fact that the variations are fast, Sfj,{t) 
being renewed after each period, while the pendulum has 
a long memory ^ e“*, i.e. many periods. So there is a 
considerable cancellation of the random component of p,. 

Slow variability however (time scale much larger than 
the period), results in little cancellation, so that the pen¬ 
dulum is quite sensitive to this type of variability. Basi¬ 
cally, we rediscover here what horologists have known for 
a long time: accurate time keeping is all about control¬ 
ling and eliminating slow, low-frequency variations in the 
parameters. Fast variations are much less detrimental. 
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FIG. 6: Simulation of the pendulum motion comprising 10® 
timesteps, or 10® dt/2'R ~ 4.8 x 10® nominal periods, about 
110 days for aim pendulum. Shown is the difference be¬ 
tween pendulum time tpp and laboratory time t (top), and 
between pendulum time ipp and model pendulum time ■i/'m 
(bottom), as a function of laboratory time. To reduce black¬ 
ening only one out of every 500 data points is plotted. Pa¬ 
rameters: iiQ = 1.5 X 10“®, e = 1.167245 x 10“"^, tpo = 0.03, 
dt = 0.03, dfJ-r.tn.s. ~ 0.3/j.o. The random component 5^(t) 
in the driving is switched on beyond the dotted line > 

4.8 x 10®). The simulation took 10“* s (double precision IDL). 


Returning to the main question, how well do Eqs. m 
and (flBl) model the pendulum dynamics? The bottom 
panel of Fig. [B] shows that the strong variability of pen¬ 
dulum time '0P seen in the top panel and induced by the 
driving noise, is by and large correctly modeled in tjjpp and 
removed. Only this small secular difference of ^ 10“® 
remains. We ran other simulations, and all graphs of 
V'p ~ V'm versus t looked more or less the same. This sug¬ 
gests that we see the effect of the 02{e^) term in Eq. (USD, 
a term that we cannot model out as its theoretical form 
is unknownii^d^ Perhaps we see the second-order cor¬ 
rection — of the oscillator frequencyj^^ That has the 
right order of magnitude and sign since —~ —10“®. 

Our conclusion is that also in this case the asymptotic 
theory appears to work well. Equations m and m 
deliver a relative rate precision of ~ e or better. This 
suffices for a physical pendulum as that attains a preci¬ 
sion of ^ £ anyway (Sec. lIV^ . In numerical experiments 
where one may switch on and off perturbations at liberty, 


there may be cases that require computation and inclu¬ 
sion of the second order terms Oi and O 2 ■ 


VII. CONCLUDING REMARKS 

The purpose of this study is to demonstrate that the 
long-term dynamics of a pendulum as typically used in 
clocks can be simulated efficiently with two first order 
ODEs, Eq. (fTSl) for the amplitude A and Eq. (ITBl) for the 
phase 'll) of the pendulum. By using the rate "0 instead of 
its inverse (the period) and by averaging over fast time 
scales, we obtain simple model equations m and (HU) 
that show that the effects of driving, friction and circu¬ 
lar error are additive. In this way we could deduce the 
expressions for these effects in a unified framework, with¬ 
out any reference to a phase diagram (Ref. [H, Ch. VII). 
In particular, the derivation of the escapement error in 
Eq. (ITOl) runs smoothly and is free of ad hoc assumptions. 

We employ the method of Bogoliubov and Mitropol- 
skji4ii5, -vvbich involves averaging over a pendulum pe¬ 
riod. Time scales smaller than a period can no longer be 
resolved, but these are not important for the time keeping 
properties of the pendulum. The pleasant consequence 
is that the model equations are rather insensitive to fine 
structure in the driving torque m(ip), which in turn al¬ 
lows us to use a simple model escapement. It is no longer 
necessary to resolve individual pendulum swings in small 
timesteps - we may now make timesteps of one pendu¬ 
lum period, resulting in a considerable time saving, and 
the possibility to make long simulations with all effects 
acting simultaneously. 

We studied two examples with the equation of motion 
o and compared the results with those of the model 
equations (HB and (HU): a sudden change in the driving 
torque and a permanently noisy driving. Satisfactory 
agreement was found. More complete models need to 
consider other disturbing factors; we only mention here 
that slow length changes Si{t) of the pendulum produce 
an extra term —51121^ on the right of Eqs. (ITUl) and (fTBl) . 
while Eqs. and m for the amplitude remain un¬ 
changed (io is the nominal pendulum length). If nec¬ 
essary new model equations may be derived for special 
escapements, see e.g. Eqs. (IC4l) and (IC5I) . 

Now that the validation has been done, the road is 
open to fast simulations with the model equations m 
and (jl6l) . For example, the simulation of Fig.|B]comprises 
~ no days for a 1 m pendulum. It took 10^ s with the 
equation of motion o, but only ~ 50 s with the model 
equations m and dUD. Note that the reliability of such 
long-duration simulations is limited by the unknown in¬ 
fluence of the second order terms. Evaluation of these 
second order terms seems to be an interesting topic now. 

It is evident from Eqs. m and (US that any fine struc¬ 
ture of the type '^n^i(flnCOsnil) -I- bnSinn'ip) may be 
added to m(ip) without affecting the time keeping. This 
opens the door to a wide variety of equivalent shapes 
for the driving torque m. But if the shape of m is so 









unimportant, why did clockmakers bother to refine es¬ 
capements? The answer must be that when parameters 
change and adapt to the environment, it is very diffi¬ 
cult not to change the two quantities j> m cos ^ dip and 
j> m sin Ip dip that are crucial for the time keeping. The 
best strategy seems to be to reduce e as much as possible. 

A pendulum is the archetype harmonic oscillator, 
an omnipresent system in physics with widely differ¬ 
ent appearancesi^ Laboratory experiments with a driven, 
damped oscillator are therefore of permanent educational 
value. These need not be restricted to a pendulum. One 
might try for example a quartz crystal or any kind of 
electronic oscillator. But whatever the nature of the os¬ 
cillator, the unifying concept of a driven, damped, rest¬ 
lessly swinging pendulum is likely to keep physicists and 
horologists spellbound forever. 


d(Asin'!/))/df = A sin')/) -|- Xipcosip. Solve for ip and in¬ 
sert A = mcosip to obtain ip = —{m/X)smip. This is 
the core of Eq. m- Add all impulses nii composing m 
in one period: S'ip = J^^iPi — “(1/'^) X)= 
— (1/A) / msinipdt (A is constant over a period). De¬ 
note the mean time derivative of ip again by ip, then ip ~ 
d'ip/2TT ~ — (1/A) J m sini/) df/27r ~ — (1/A) j> m sin ip dip 
(since dip/dt = 1-|-small terms, so dt ~ di/). Finally, add 
the nominal rate ip = 1, and the circular error —A^/16. 

It is hoped that these informal arguments help to 
clarify those of Ref. [l|, Ch. 7 and Woodwardj^i 


Appendix C: A single-impulse escapement 
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Appendix A: The circular error 


The circular error term follows from Eq. (|S]): 


circ. eiTOi — j){X ^ sin(Asin'!/) — sini/} sin^ d')/ 


Here, Ji is a Bessel function of the first kind<^ The es¬ 
sential steps in the computation above are (a) to use for 
sin(Asin'0) relation 9.1.43 of Ref. IXl and then (b) the ex¬ 
pansion Ji(A) = x(l — jVl -) with X = A/2, 

cf. relation 9.1.10 of Ref. 

We present Eq. (EU here because it differs from the 
usual expansions (e.g. Ref. m, p. 51 and Ref.0, Eq. (8). 
This is because Eq. (EH) specifies the effect of circular 
error on the rate, not on the period. 


The model escapement, Fig.jH that we have used sofar 
interferes continuously with the pendulum motion. Here 
we compare it with an escapement of opposite type, that 
delivers in each period a single impulse at a fixed pendu¬ 
lum angle ^p\. 


= eS{ip - ipi) . (Cl) 


The function S{x) is basically a sharp peak at x = 0. 
Width and height are inversely proportional such that 
J S{x — xq) dx = 1, provided Xq G I, the integration in¬ 
terval. The width of S{x) is supposed to be much smaller 
than all fine structure in the problem at hand, whence 
/ f(x)5{x — xq) dx ~ fixo) if xq G I, otherwise zero. 

The energy transferred to the pendulum in one period 
is J m dip = e > 0. It follows that dip > 0 ad ipi, i.e. ip is 
increasing. Since ipi = Asini/i according to Eq. the 
phase ipi at the impulse is in the 1st or 4th quadrant. 
We compute a and b given by Eqs. cm and m- 


b = (p msinip dip = — / S {ip — ipi) sin ip dip . (C2) 

J 27r 


Here we hit a technicality: to be able to apply the S- 
function recipe we transform the integration over ip to 
one over ip. From ip = Asin^/ and A ~ constant in a 
period, we infer dip = {dip/dip)~^dip = (Acos^)“^d(/3, so 

C 6 

^^XtAJ t,anip 6{ip - ipi)dip =-^^tanipi . (C3) 


Appendix B: Back-of-the-envelope derivation of 
equations dSD and (nni) 


A similar computation yields a = ej (27rA). To obtain the 
new model equations for A and ip we repeat the calcula¬ 
tion of Sec. IHIBl resulting in: 


The power fed into the pendulum is torque x angular 
speed = mip = mX cos ip (use Eq. (ED for ip; gravity does 
on average no work). But this is also the time derivative 
of the total energy, so d(^A^)/dt = mXcosip or A = 
mcosip. Eq. ED follows by adding friction and averaging. 

Suppose m is a single impulse. At that instant ip 
gets a boost while ip = A sin is constant, so 0 = 


A = 

Ip — 1 = 


A2 e 
16 27rA2 


tanipi . 


(C4) 

(C5) 


These equations may be used to study the performance 
of a single-impulse escapement. We briefly look at the 
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case of stationary operation. The amplitude is then A = 
and the rate is: 

- 1 = -etanV'i ■ (C6) 

Al-ne 

The last term is the escapement error, which now either 
increases the rate {ipi < 0, impulse before a zero), or 
reduces it {ipi > 0, impulse after a zero), cf. Ref. [ll, p. 
130. Note that ipi = Asin'i/ji implies tan'i/ji = — 

and just as in Eq. (fTCl) . we see the abearance 
of the square root functions in Ch. VII of Ref. [l|. 

Appendix D: Numerical issues 

We rewrite Eq. m as a first order ODE in the vari¬ 
able {(fi, (fi) which we solve with the IDL routine RK4. 


To model the escapement the code finds out, at each 
timestep, in which quadrant ip lies (from the signs of ip 
and ip). From the values of ip and ipo the code then de¬ 
cides between m = +pL and m = —pt, according to Fig. [5] 


The simulation delivers pendulum angles ipk on an 
equidistant grid tfc, stepsize dt. To measure the rate 
with Eq. (HH) we need the —h zero crossing times of 
ip with a precision much better than ^ e. This follows 
from Eqs. (ED and (1^^ . It is therefore necessary to 
do the whole (IDL) simulation in double precision. The 
timestep k enclosing a zero obeys ipk-i < 0 and ipk > 0. 
An approximate zero crossing time is found by linear in¬ 
terpolation: 4ero = tk -1 + dt ■ ipk-ll{ipk-l - p>k)- This 
gave an acceptable accuracy, though we do observe some 
numerical noise, for example in Fig. [SJ top panel. 
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